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ABSTRACT 

We develop a diagrammatic technique to represent the multi-point cumulative prob- 
ability density function (CPDF) of mass fluctuations in terms of the statistical prop- 
erties of individual collapsed objects and relate this to other statistical descriptors 
such as cumulants, cumulant correlators and factorial moments. We use this approach 
to establish key scaling relations describing various measurable statistical quantities 
if clustering follows a simple general scaling ansatz, as expected in hierarchical mod- 
els. We test these detailed predictions against high-resolution numerical simulations. 
We show that, when appropriate variables are used, the count probability distribution 
function (CPDF) shows clear scaling properties in the non-linear regime. We also show 
that analytic predictions made using the scaling model for the behaviour of the void 
probability function (VPF) also match the simulations very well. We generalise the 
results for the CPDF to the two-point (bivariate) count probability distribution func- 
tion (2CPDF), and show that its behaviour in the simulations is also well described 
by the theoretical model, as is the bivariate void probability function (2 VPF). We 
explore the behaviour of the bias associated with collapsed objects in limit of large 
separations, finding that it depends only on the intrinsic scaling parameter associated 
with collapsed objects and that the bias for two different objects can be expressed as 
a product of the individual biases of the objects. Having thus established the validity 
of the scaling ansatz in various different contexts, we use its consequences to develop 
a novel technique for correcting finite-volume effects in the estimation of multi-point 
statistical quantities from observational data. 
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1 INTRODUCTION 

The process of gravitational instability leads to the growth and collapse of initially small inhomogeneities in a self-gravitating 
medium. This is the basic idea behind theories for the formation of galaxies and large-scale structures in the Universe, as the 
small initial fluctuations seed the formation of collapsed objects that progressively cluster and merge into objects of larger 
size, such as galaxies and clusters. Much progress has been made during the last two decades using numerical simulations 
to follow the evolution of cosmological density fluctuations up to and beyond their collapse into bound structures, but it 
is nevertheless important to develop an analytical framework for describing this phenomenon and understanding what the 
simulation results. 

Much of the previous analytic work in this direction has been based on an extended version of the Press-Schechter (1974) 
model. Such methods have been used to study the distribution of formation epochs (Bond et al. 1991; hacey & Cole 1993; 
see White 1993 for a complete review), merger rates and survival times of haloes and also how the distribution of collapsed 
objects relates to the underlying mass distribution (Mo & White 1996; Mo et al. 1997). Tests of these predictions have shown 
good agreement with theoretical model predictions. Other recent methods which try to explain statistics of collapsed dark 
objects depend on extension of Zel'dovich approximation (bee & Shandarin 1997, 1998; Catelan et al. 1998; Porciani et al. 
1998). 
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A parallel approach has been to assume that statistical properties of the distribution of collapsed objects - particularly 
the many-body correlation functions - follow a hierarchical scaling pattern (Balian & Schaeffer 1989a, b). This kind of scaling 
can be used to predict, for example, the mass function of collapsed objects (Valageas & Schaeffer 1997) and their correlation 
properties (Bernardeau & Schaeffer 1992; Munshi et al. 1998b). These results are especially interesting because they relate 
directly to the statistics of the underlying mass distribution and hence to the dynamical origin of the distribution governed 
by the BBGKY equations. 

Predictions of the behaviour of one-point statistical quantities such as the count probability distribution function (CPDF) 
or void probability distribution function (VPF) for the underlying mass distribution in the highly non-linear regime of strong 
clustering were derived using this kind of scaling ansatz by Balian & Schaeffer (1988, 1989) . These results were later extended to 
the quasi-linear regime where a perturbative series expansion can be used (Beranardeau 1992, 1994a, 1994b, 1995). It emerged 
from these and other studies that a scaling hierarchy naturally developed by the action of gravitational clustering on Gaussian 
(random-phase) initial perturbations. Predictions for these one-point quantities have been tested with reasonable success 
against observations using various galaxy catalogs (Maurogordato & Lachieze-Rey 1987; Blanchard et al. 1990; Mourogordato 
et al. 1992; Bouchet et al. 1993) and numerical simulations (Bouchet & Hernquist 1992; Lucchin et. al 1994; Colombi et al. 
1992, 1994, 1995, 1996; Munshi et al. 1998c), although Lucchin et al. (1994) and Colombi et al. (1996) found some deviation 
from scaling. 

It is natural now to ask whether these properties of one-point statistical quantities can be generalised to two-point 
statistics, or even beyond. The two-point count probability distribution function (2CPDF) and two-point void probability 
function (2VPF), both defined in detail below, are particularly interesting for a number of different reasons. The first is that 
the 2VPF can be shown to be a generating function for the cumulant correlators (CCs) in a way which is very similar to 
the VPF, which acts as a generating function for cumulants. Recently it was shown that CCs can be reliably measured using 
galaxy catalogs (Szapudi et al. 1992; Meiksin et al. 1992; Szapudi et al 1995; Szapudi & Szalay 1997). They were also measured 
in numerical simulations (Munshi & Melott 1998) to test extended perturbation theory (EPT) as proposed by Colombi et al. 
(1996). Since it becomes increasingly difficult to measure these quantities with increasing order, it is sensible instead to study 
their generating function (2VPF) which carries information to all orders. The second reason is that these two-point objects 
are useful in the problem of estimating errors associated with the determination of one-point quantities such as the CPDF, 
VPF and Sn parameters (Szapudi & Colombi 1995). Although analytical results related to two-point quantities are frequently 
used to estimate and correct such errors, no systematic studies have yet been done to test these predictions directly against 
numerical simulations, partly because of the difficulty of performing numerical simulations with the very large dynamic range 
necessary to study these quantities. Finally, the 2CPDF carries information as to how overdense regions are biased tracers of 
the underlying mass distribution, so a computation of the bias implied by the scaling ansatz can provide an understanding of 
the origin of bias that can be compared with results from other approaches. 

Before proceeding, however, we stress that one should not get too carried away by the success of this kind of scaling model. 
The hierarchical ansatz is, at best, a simplifying hypothesis which, admittedly seems to work very well for "realistic models" 
of structure formation. But no firm connection has yet been established with the microscopic physics of collisionless clustering 
described by BBGKY equations (Peebles 1980). Most efforts in solving the BBGKY equations have focused on particular 
closure schemes (Davis & Peebles 1977), the general separability of position correlations and momentum correlation in phase- 
space (Fry 1984; Hamilton 1988) and its stability properties (Ruamsuwan & Fry 1992; Yano & Gouda 1998). Although these 
precious efforts in solving the highly complicated nonlinear integro-differential BBGKY equations provide us with valuable 
insights, they also help underline the complications in arriving in any general solution to the problem. Our inability to solve 
these equations in the highly non-linear regime leaves us with undetermined parameters which can only be computed from 
numerical simulations, thus complicating the task of testing the scaling ansatz itself. 

This paper is structured in the following way. Section 2 is devoted to a theoretical discussion of multi-point cumulants, 
multi-point factorial moments, multi-point factorial correlators and their various generating functions. We also develop a 
general approach to the treatment of finite-volume effects using these quantities. Section 3 is devoted specifically to one-point 
and two-point statistics. We summarize the main analytical results of Balian & Schaeffer (1989) and Bernardeau & Schaeffer 
(1992) and test them numerically. Section 4 contains a detailed description of our numerical simulations and to details of data 
analysis techniques; we also present the results of our analysis in this section. We discuss the main results and relate it to 
other studies in Section 5. 



2 VOIDS & PROBABILITIES FROM CUMULANT CORRELATORS 

The void probability function (VPF) is the probability that a randomly-placed sphere of some volume V contains no galaxies. 
This can be written 

P.(0)=exp(-^)), (1) 

in which the mean number of galaxies per cell is N and the volume-average of the two-point correlation function over the cell 

is 

? 2 =V72 i i Un,r2)dV 1 dV2. (2) 
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The quantity N c — N£ 2 - The VPF is a generating function of P(N), the count probability distribution function (CPDF) for 
cells of the same size (Balian & Schaeffer 1989). Let us assume Using the generating function relationship between the CPDF 
and VPF we can write 



which can also be written in the following integral form (Munshi et al. 1998c) 



(3) 



The integral is to be evaluated along a contour around A = 0. We also define the function Vl/(t) = —<j>(—t) which can be 
directly related to the function n(^), the continuous analogue of P(N), 

exp (^)=r* /exp (^) n(,/) - (5) 

Using the definition of P(N) now one can easily show that 

P(N)^J\,^t^n { u). (6) 

This is the form of distribution obtained when II(i/) is Poisson sampled, because it is the convolution of the function II(^) 
with a Poisson distribution describing the shot noise effects. This demonstrates that II(i>) can be viewed as the continuous 
limit of P(N) in the limit of large number densities. This can be seen by change of variable A = 1 + t/N c in equation (3) and 
then taking the limit N c — > oo and N — > oo with the ratio N/N c remaining finite; this also gives P(N) = H(N). 

The factorial moments of P(N) can now be related to the moments of il(z^) with the aid of the definition of the p-ih 
factorial moment: 



OO />OC 

jv-n -A> 



l{v)dv. (7) 

One can similarly define the normalized factorial moments of P(N) by 

E » = pjk £ N ( N - V"< N -P+ V P W = pjk (N ^- ^ 

One can define a generating function E(t) for the S p parameters using 

oo 

S(t) = ^S p ^ (9) 

i 

and we have the following expression connecting these two generating functions ^(t) and £(£) (Munshi et al. 1998c), 

l + f)=exp ( m, (10 ) 

?2 ?2 

which can also be written as 

*(t) = £ 2 ln(l + H2). (11) 

?2 

It is possible to extend the method of generating functions to multi-point factorial moments. If we denote the two-point 
CPDF P(Ni,N 2 ), the joint probability of finding Ni particles in the first cell and N 2 particles in a second cell, and assume 
locally Poisson sampling distribution then we can relate P(N\, N 2 ) with its smooth counterpart tt(ui,v 2 ) by: 

/>oo />oo / \ N-i 1 \ No 

P(N U N 2 ) = duj „ 2 exp(-J> 1 exp(-^ 2 



Nil N 2 \ 



so that the two-point factorial moment becomes 

00 poo poo 

{(TViA^M = J2 ^1(^1 - 1)-- ■ (M -P+1)N 2 (N 2 - 1)... (JV a -q+l)P(N u N 2 ) = / v\dv v / v\dv 2 U(v u v 2 ) (13) 



AT=0 



JO 



We define the normalised two-point factorial moment by following equation. 

F2 c'2 



<2) 



s # = , ,1%+, y)Wi(Wi - 1) . . . (JVi - p + l)^(iV 2 - 1) . . . (JV 2 - q + 1)P{N U N 2 ) = , ^ rP+q (NiN 2 ) pq (14) 

p\q\Nc — piq'.Nc 

It is to be noted that, in general, owing to statistical isotropy and homogeneity S pg should depend only on the separation of 
two cells. We define the generating function for E pq by E*- 2 ', i.e. 
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oo 

E (2) = ]T E M ^ (15) 

and similarly the generating function for the 2CCs by 

£2 Z^ p\ q \ (,52( a; ))(P+9) • l j 

Using similar technique as used for the case of one-point moments we can relate E' 2 '(£i,£ 2 ) with \I/' 2 '(ii,£ 2 ): 
E(ti) E(t 2 ) S^fo.ta) _ /*(ti) *(t 2 ) gWfti .ta)' 



x + -v^, + -^ + - ^.^ = cxp i^n + i^i + - ^»*> (17) 



A Taylor expansion of the above equation provides us with a set of relations which can be used to measure 9pJ once E p( / is 
computed: 

*£? = ("l + Ei 2 ?) (18) 

*!? = i(l-EW+6(E^-E 2 ) (19) 



6 



2 

*<? = i(-3 + 4E<? - E< 2)2 - 2£ 2 (E< 2) - 2E 2 + E 2 2) ) - 2& 2 (E 2 - E 2 2) )) (20) 

S2 

*£> = i(-l + Ei?-6(sg ) -2E a )-6Ei?E 3 + g(E«-E 3 )) (21) 

S2 

*£> = l(2 + Ef 1 )2 -E 11 (3 + 6(Ei 2) -2E 2 )) (22) 

S2 

+£ 2 (2E< 2) - 4E 2 + E 2 2) )) - |f (E< 2) - 2E^ + S 2 S 2 2) + S 22 - S 3 ) + |f (E 23 - E 2 E 3 )) 

In our above analysis we have used two different cells of the same size, but it is now obvious that the corresponding result for 
two cells of different size will read 

Efti) S(fa) S (2) (*i,frO = ( 9(h) tt(* 2 ) *< a >(ti ,fe) 

tx £W ^z^" I £z £2/ tit' 

S 2 S 2 S 2 S 2 \ %2 S 2 S 2 S 2 

These results can be extended reasonably straightforwardly to three-point quantities 



1 + t^ll + ^11 + - ^11 = exp ll^J. + 1^1 + ~ v-i-^ (23) 



E(ti) , E(t 2 ) , E(t s ) , EW(ti,fa) , £< 2 >(ti,t 3 ) , E< 2 >(i 2 ,i 3 ) , E (3) (ti,t 2 ,*s) 
where 



?2 ^2 ^2 ^ 2 ^ 2 ^ 2 ?2 



*> = -?— + -?— + -j— + ~ 2 + T2 + T2 + ?S ( 25 ) 



tt(ti) £(fa) £(fa) ^ (2) (ti,t 2 ) ^ (2) (fa,t 3 ) ^ (2) (ti,fe) tt (3) (ti, t 2 ,t 3 ) 
la 6 6 CI g g CI 

and 

£3 Z^ p!g! r ! £<P+9+r) *■ > 

The quantity E' 3 ' is similarly the three-point generalisation of the two-point quantity E^ 2 '. Expanding this relation one can 
express three-point moments in terms of factorial moments of the three-point CPDF. 

*g\ = (2-3E( 2) +E< 3) 1 ) (27) 

*ii2 = r (~ 3 - E n 2 " E m + l» E ua ~ ^i? + 2 ^ Sa + E ii ) ( 5 ~ & Ea )) (28) 

S2 

*8s = i (4 + 2E 2 ! + Eg\ + 4&E< 2) - 2gE<|> - 6&E 2 - f 3 Eg\ E 



f 2 

^2 



2 



-6Sg a + 2|f E 3 - E< 2) (7 + 2C 2 (E( 2) - 3E 2 ) + |f E 3 ) + gEn 3 ) 
It is possible to generalise these equations for an arbitrary number of points: 

HEf + E^^ E S(3)( y' tfc) +... + S(P)( |- U) -ex PXp , (29) 

pairs " triplets 



where 

S2 . . , . 

pairs triplets 



Xp -2-^r + ^ i + ^_ I ••■ + i • (30) 
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It is important to note that two-point quantities depend only on the separation of two cells, but the multiple-point correlators 
and their generating functions depend on the geometrical configuration of different cells. We have used the following definition 
for multiple-point cumulant correlators: 

*W(ti,fa,...t fc ) ^ ( t p 1 t%...tl \ (8"(xi)S^x 2 )8 r (x 3 )...S s (x k )) 

>2 



£* £-~t yplqlrl ...s\ J c 



Inverting this relation, we can express ^^ k \t\,t 2 , . . . , t k ) in terms of the generating function for factorial moments, which can 
be measured directly from simulations: 

* (fc) (ti,t 2 ,...,i fe ) = lnM k (ti,t 2 ,...,t k )~ Y lnM fc _i(ti,t 2l ...,tfc-i)+ Y M k - 2 (ti,t 2 , . . . , t k - 2 ) + (32) 

(k — l)tuplcts (fc — 2)tuplcts 

... + (_l)( fc - 2 )^lnM 2 (ti,t i ) + (-l) ( *- 1) Y^ In M(U) 

pairs singlets 

where we have defined 

M(ti) = 1 + ^ (33) 

?2 



S(fa) S(fa) T.^{ti ,t 2 ) 

6 6 CI 



Mafr.fa) = l + ±^Z + ^ + - ^ . (34) 



It is also possible to work with factorial correlators instead of the factorial moments defined above by extending the two-point 
factorial correlators introduced by Szapudi & Szalay (1995) to multiple points. If we denote the k th order factorial moments in 
different cells by (Ni) k then we can write the generating function of the two-point factorial correlators W k i by W(ti, t 2 ), the 
generating function of the three-point factorial correlator W k i m by W(ti, t 2 ,tz), and so on. In general W(ti, . . . ,t 3 ) denotes 
the generating function for s-point factorial correlators W k ... r : 



WWfr) = ! + £-(£) {(N l)k ) (35) 

fc=l 

W {2 \ tl ,t 2 ) = Y M (|) (|) (((Ni)* - m) k ))((N 2 )i -{(N 2 )i))) (36) 

00 , . 

W^iUMM) = Y fcilLld) (I) (|) m (((^) fc -((A r i) fc »((A r 2) i -<(iV 2 ) i ))((iV 2 ) m -((iV 2 ) m ))) (37) 

fc=l,l=l,m=l 

wM(t u ...,t.) = Y fcT^nd) ■■■(^) r (((Ni)k-((N 1 ) k ))...((N 3 ) r -{(N s ) r ))) (38) 

fc = l,...,r=l 

Now it is possible to link the basis function W^-"' for computing multi-point cumulant correlators with the other set of basis 
functions E' n ', which carries equivalent information. 

= i+m 09) 

£2 

= lE (2) (t 1 ,t 2 )-ls(t 1 )ls(t 2 ) (40) 

?2 ?2 42 



w ( 


vf 2 / 


IV (2) (^. 

^?2 


6/ 


(3)/£l *2 


>1 

6/ 



1 "(3) I .... ' .... I 

!„,. , l„m,. . , 1„,. , 1 „ m ,. . , 1„,. . 1 



-3E w (ii,i 2 ,t 3 ) + 2— E(ii) — E(£ 2 ) — E(i 3 

?2 ?2 £2 £2 



E(ti)^E«(fa,t 3 ) - -S(t 2 )- f E (2) (t 1 ,t 3 ) - E(i 3 )- 2 E (2) (t 1 ,t 2 ) (41) 

?2 4 2 ?2 £ 2 £ 2 t, 2 

Relating multi-point factorial moments with generating functions of multiple-point cumulant correlators (defined earlier) we 
can write: 

exp(xi) = ex p(^)=^(j|); ( 42 ) 

«p(»» - « P (iM + iMli^M)^^(|,|) + ,,(|),,(|) 

«p(„) = ff«.(|,|,|) (44) 

+H " 2 4f Ml) +H " 2 '(|<f Ml) + M|<f Ml) M|M|M|) <«> 
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Generalizing these results, and connecting multi-point factorial correlators with multi-point void probability functions we 
obtain the result equivalent to eq(E9l). 

»""(| |)+ e """(D^'-'d ^ )+£»'(sWtK- 2 '(! i?) <*» 

singlets pairs 



+ ...+ £ lV (|)^ W (i |f )M,<»(|.|) + »-(|)..,,(| 



(s — 2)tuplcts 



>p ^ (3) (^,M ^ (p) (^i^2,...i fc ) \ 

^ el " e, r 



pairs triplets 

It is to be noted that multi-point cumulant correlators can also be defined in such a way that they incorporate all lower order 
correlations. For example, if we allow one or more of the powers p, q, r or s to vanish the correlator becomes independent of 
one or more spatial co-ordinates. This will allows a more compact way to express factorial moments in terms of multi-point 
cumulant correlators (Szapudi & Szalay 1997) which is equivalent to the result obtained using factorial moments. However, 
it is useful to decompose multi-point factorial moments and multi-point factorial correlators in terms of irreducible one-point 
and two-point quantities from a computational point of view. As we shall see, this also allows us to relate multi-point cumulant 
correlators directly with the tree amplitudes involved in the hierarchical ansatz. 

The discussion so far is completely general. While Szapudi & Szalay (1993) used the hierarchical approximation of 
Qn+m — QnQm and analogous approximations for higher orders, we will be concentrating on the formalism developed 
by Bernardeau & Schaeffer (1992) in which amplitudes associated with hierarchical trees are assumed to be the product of 
tree vertices. Following paper I, we can now relate ^' n '(ti, . . . ,t n ) with generating functions ix n {t) of vertices appearing in 
tree-representation of multi-point cumulant correlators which depends on the hierarchical ansatz. These vertices depend on 
the tree-level amplitudes v n in representation of correlation hierarchy of matter correlation function (see paper I for complete 
discussion): 

*(*) = E^r (47) 

71=1 

n = l 

M3(f) = J2^f; (49) 

n=l 

etc. The quantities /i n determine the amplitude associated with a vertex that has n external legs. Using the vertex generating 
function we can express & n '(ti, . . . ,t n ) as: 

* (2) (ii,i 2 ) = jUi(-ti)U/ui(-ti) (50) 

* (3) (ti,t2,t3) = Hl(-ti)£ a bH2(-ti)£aclJ<l(—ti) + • • • (cyclic permutations) (51) 
* (3) (ti,i2,i3,i4) = Mi( - *i)£ab/43(— £l)£a.eMi(— *i)£ad/^i(— tl) +••• (cyclic permutations) 

+Hi{-ti)£,abLJ.2(-ti)£;acH2(-t 1 )£adLJ.i(-ti) + . . . (cyclic permutations) (52) 

In deriving above relations we have considered only the dominant contribution and this approximation is valid only when the 
density variance in each cell is much higher than the correlation between cells. Szapudi & Szalay (1997) proposed the use 
of the two-point cumulant correlator to separate amplitudes associated with different tree vertices, but since the number of 
tree topologies increases exponentially with order of the diagram (Fry 1984), the two-point cumulant correlator can only be 
useful in up to 4 th order diagrams. This deficiency can only be cured by moving to multi-point cumulants and our formalism 
developed here will be able to determine v n parameters for arbitrary n. However, it is obvious that multi-point cumulants 
are in general dependent on the geometric configuration of the points, and their determination becomes more complicated as 
the order increases. We also stress the point that our formalism can in general be used in principle to determine all the tree 
amplitudes without making any assumptions at all. This is very useful as this can actually test the validity of all hierarchical 
models of gravitational clustering in the highly nonlinear regime. 

The multi-point CPDF (MCPDF) can be decomposed in terms of several one-point quantities. Expanding the multi-point 
void probability function, which acts as a generating function for MCPDF, Bernardeau & Schaeffer(1992) showed that it is 
possible to express the MCPDF in terms of a series expansion in £ij/%2 « 1, where £y is the two-point correlation function 
between different overdense cells and £2 is the variance within one cell as defined above. Such a decomposition is expected to 
work for overdense cells when they are separated by moderately large distance. To first order in such a decomposition it was 
found that 

P(N 1 ,N 2 ) = P{Ni)P(N 2 ) +P(7Vi)b(A f i)Ci2(ri2)P(iV2)fe(A r 2), (53) 

where b(Ni) is the bias function for overdense cells which describes how over-dense cells are correlated compared with 
background mass. Normalisation constrains on P(Ni,N2) can be translated in to constrains for b(N). 

© 1996 RAS, MNRAS 000, 000,000 
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P(N) 

P(N)b(N) 



P(N)V 2 (N) 

P(N)V(N) 
3 



MOMENTS 
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^P 
v v 3 



P(N)V(N) 
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P(N,. 

l 



N) 



Figure 1. Basic elements in the representation of multi-point joint CPDF and related factorial moments. We represent P(N) and one 
associated factorial moment E p by a circle with no legs; P(N)b(N) and associated moments Ep 1 are represented by circle with one leg 
and, similarly, a circle with n legs represents P(N)u„(N) and its associated factorial moments £p", where p represents the order of the 
relevant factorial moments, and n the order of the relevant vertex. 



J2KN)P(N) = 0; 
V^AT6(iV)P(iV) = N. 



(54) 



With such a decomposition we can express two-point quantities, such as Ep 9 , in terms of one-point factorial moments but 
with the effect of bias b(N) taken into account: 



2-*pq — Z-ipZ-iq T" 2-*pS12^q, 

where we have defined a new biased one point factorial moment Ep by, 

6 



E ? = 2 7m? N{N ~ 1 )---( n -p+ mN)P(N) 



pW, 



(55) 



(56) 



A similar decomposition can be performed for the 3-point CPDF and, if we write the expression up to second order in (,ijS, 2 , 
we get 

P(N U N 2 ,N 3 ) = P(N 1 )P{N 2 )P(N 3 ) (57) 

+P(N 1 )P(N 2 )b(N 2 )^ 3 P{N 3 )b(N 3 ) + P{N 2 )P(N 3 )b{N 3 )^iP(N 1 )b(N 1 ) + P{N 3 )P{N 1 )b{N 1 )^ 2 P{N 2 )b(N 2 ) 

+^(Ar 1 )P(Ar 1 )ei 2 K^2)P(A r 2)Ci3fc(A r 3)P(A r 3) + ^(A r 2)P(A r 2)63&(A r 3)P(iV3)C2iK^i) J P(^i) 
+i-2(N 3 )P{N 3 )^i3b(N 1 )P(N 1 )^ 2 b{N 2 )P{N 2 ). 

A new function appears, together with its associated constraints: 



^is 2 (N)P(N) = 
^2Nv 2 (N)P(N) = v 2 N; 



(58) 



v 2 is the second-order vertex associated with a tree-development of the correlation hierarchy in the highly non-linear regime 
We define the associated one-point factorial moment Ep, which takes into account of corrections second order in £ij£ 2 : 

6 



E p 2 = S plk N{N -!)••• (^ - P + l)^2(iVi)P(iVi). 
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Figure 2. The decomposition of the MCPDF is represented diagrammatically. Every straight line joining two spheres represents the 
correlation between two cells. Higher-order terms contain more straight lines joining a given set of cells and hence can be ignored in a 
leading-order approximation. Only tree terms are shown. All possible permutations of labels are to be considered, but in this figure we 
show only one diagram for each topology. In second order in £;j /£; , it is possible only to have one type of topology, in third order we can 
have a two different topologies represented by "snake" and "star" diagrams. In even higher order it is possible to have a hybrid of these 
two different types. It is interesting to note that the correlation hierarchy for matter induces exactly similar hierarchy for MCPDF. 

With this higher order contributions are taken into account E pgr can be expressed in terms of one point quantities such as 

2-1 p . 2-1q cIIICI 2-Iq , 



EpSgS r + SpS g ^23S r + £ g £j,£l3£ r + 2 r £j,£l2£ q + £g 2 £23£p£l3S r + S^ 2 £3i£j,£l22 g + Sp 2 ^3iS p ^32S g 



(60) 



It was shown in paper I that such a decomposition is possible to arbitrary order, as long as we neglect the contributions 
from loop terms. The tree structure of MPCDF and associated multi-point factorial moments is induced by the inherent 
tree structure of correlation functions and is very similar in nature. It is interesting to note that, in this particular case, 
contributions from different orders are completely separated and also the contribution to joint moments from different spatial 
locations can be expressed as different independent local terms depending only on one spatial co-ordinate. This also makes 
volume corrections easier to perform, as we shall see in §4. 

3 SCALING IN ONE-POINT AND TWO-POINT QUANTITIES 

Scaling properties that arise because of the inherent tree structure of correlation functions in highly non-linear regime were 
studied by Balian & Schaeffer (1989). The principal assumption in this and related studies is that the vertex amplitudes 
which appear in a tree development of higher-order correlation functions should be constants. Every higher-order correlation 
function is build by taking suitable products of two-point correlation functions and all such N different tree's carry different 
amplitudes: 



&v(ri 



. rjyl 



_ __, (jv_1) 

2 , Tjv - a 2^ II kri,Tj)- 

a, N — trees labclings edges 

which in most general case can be expressed as 
fiv(Ari, . . . Arjv) = A~' T(iV_1) £jv(ri, . . .r N ). 



(61) 



(62) 



Clearly no loop terms (such as those that appear in the Kirkwood scaling relation) are considered in such an analysis. There 
is some support from observations for such an assumption in highly nonlinear regime. On the other hand, in the quasi-linear 
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regime and in the limit of vanishing variance such a hierarchy occurs naturally but the hierarchical coefficients in general 
depend on different shape parameters appearing in higher-order correlations. Recent studies of the bispectrum, however, show 
that nonlinearity tends to smooth out such shape dependence and in the highly non-linear regime the bispectrum becomes 
completely independent of shape. Such a result of course is in agreement with the hierarchical ansatz where all higher-order 
parameters are assumed to become shape independent (Scoccimarro et al. f998). 

Any assumption regarding the entire hierarchy of many-body correlation functions has very strong consequences and 
can be used to derive properties of one-point and multi-point statistics, ft is important to note that since the exact values 
of hierarchical coefficients are left unspecified in such an ansatz, it is only possible to talk about generic scaling properties 
induced by them. 

The void probability function (VPF), which is very simple to measure from N-body data or galaxy catalogs carries 
information on volume averages of all orders of correlation functions. It has been shown that VPF -fV(O) which denotes the 
probability of finding a cell of volume V empty of any particle can be expressed as an unique function of scaling variable N c , 
defined above in §2(White f979; Balian & Schaeffer f989 ). In particular we have (Balian & Schaeffer 1989) 



iV(0) = exp (-Na(N c )) = cxp (~^f^) 



= cxp ( =- ) (63) 



For large values of N c , a(N c ) becomes a power-law with power law index —u) that depends on initial conditions 

a{N c )= aN~ w (64) 

(Balian Schaeffer (1989). Note that the VPF depends by construction on the discrete nature of sampling and it is not possible 
to define a VPF for a continuous matter distribution. In effect, however, the VPF can be used to learn about the Sn parameters 
of the underlying field. 

Scaling in VPF induces a similar scaling in the CPDF, and the evolution of the CPDF through different length-scales 
and different epochs can be described using a unique scaling function h(x) (with x = N/N c ). This function is related to the 
VPF through 

/ioo , 
|£*(y)exp(ys) (65) 

■ZOO 

(Balian & Schaeffer 1989) and Py(N) is related to h(x) by the following expression (Balian & Schaeffer 1989), 

^ (iV) = ^(f) (66) 

Using these relations it is now possible to express Pv(N) as 

Pv{N) = mTn(N:) ^ 

(Balian & Schaeffer 1989). This expression is valid in the range N v < N < N c where, physically, N v is the typical cell 
occupancy in underdense regions and N c in overdense regions; these are related to each other by N v — AT c (a/£2) ■ It is 

interesting to note that the power-law index — u), which appears in the expression of the VPF, also appears in case of CPDF. 
For highly overdense cells a different asymptote can be computed by noting the that <j>(y) exhibits a singularity for small 
negative values of y = y s 

4>(v) = <Ps - a s T(io s ) (y-y s )~" s (68) 

It introduces an exponential cut-off for large values of JV in P(N) (Balian &: Schaeffer, 1989), 

h(x) = a s x^ s ~' exp(—\y s \x), (69) 

Previous studies have shown that an analytical fitting function for h(x) works well in both these regimes (Balian & 
Schaeffer 1989): 

The constraints satisfied by h(x) are 

© 
xh(x)dx = 1 (71) 

o 

and and 

x 2 h{x)dx = 1. (72) 

'0 

The hierarchical ansatz has also been used to predict a scaling behaviour of the multi-point CPDF and related statistics. 
The formalism developed by Bernardeau & Schaeffer (1992) shows that the two-point CPDF can be written as 
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Pv(N u N 2 ) = P(N 1 )P(N 2 ) [l + €a(ri,ra)6(JVi)6(JV 2 )] , (73) 

where we have denoted joint occupation probability of two different cells (we take sizes of these two cells to be the same) by 
Pv(Ni, N2). Such a decomposition means that the correlation function of two cells with occupancy Ni and N2 (Bernardeau 
& Schaeffer 1992) is 

6(JVi,JV 2 ) = 6(Ari)6(JV 2 )€a(ri,r 2 ). (74) 

Several characteristic features of the bias induced by these scaling properties can be noted. For one thing, it is clear from 
this derivation that the bias depends only on x and hence only on intrinsic properties of the collapsed object, at least for 
large separations. Since the bias b(N) satisfies the property ^2 N P(N)b(N) — it is clear that if the bias for some occupation 
numbers is positive, it must be negative for others. Detailed calculations show that cells with occupation number close to N v 
(the average occupation of cells in underdense regions) have negative bias. On the other hand, the correlation of two identical 
objects (i.e. cells with same occupation number) is always positive: 



b(x) = — ; (for N>N C ) (76) 

Interesting though b(x) is in itself, from a computational point of view it is much easier to compute 6(> x), the bias of 
overdense cells above certain threshold, defined as 

f °° b(x')h(x')dx' 

b(> x) = Jx ' . (77) 

f™h(x')dx> 

It can be shown in the two regimes we have discussed that 6(> x) can be related directly with b(x), i.e. b(> x) = 2b(x) 
for N < N c and 6(> x) = b(x) for N > N c . It is interesting to point out at this point that for a Gaussian density field 
b(> x) = b(x) — 1 for rare events. In the quasilinear regime, the value of b(> x) computed numerically has been shown to 
match theoretical predictions (Bernardeau 1996). However no such studies have been made in the highly nonlinear regime 
where we expect the hierarchical ansatz to be valid. 

In paper I we derived higher-order predictions of the hierarchical ansatz, but we restrict ourselves here to two-point 
quantities associated with overdense cells. The bias in this case acts as a generating function for the cumulant correlators; 
this was studied in Munshi & Melott (1998). A complete test of higher-order predictions covering the CPDF and VPF of 
overdense cells will be presented elsewhere. 



4 SIMULATIONS AND DATA ANALYSIS 

4.1 Simulations 

The simulations used here are numerical models for the gravitational clustering of collisionless particles in an expanding 
background. We study evolution of initial two-dimensional Gaussian perturbations in 57 = 1 universe. All the 2D simulations 
are done with a particle-mesh (PM) code with 4096 2 particles with an equal number of grid points. The code has at least 
twice the dynamical resolution of any other PM code with which it has been compared. The 2D simulations we use here are 
very similar to those described in detail in Beacom et. al (1991), a video of their evolution in shown in Kauffmann & Melott 
(1992), but with 64 times more particles. In this paper we analyze a subset of the simulations with featureless power- law 
initial spectra of the general form, 

P(k) oc k n for k < k c , 

= for k> k c . (78) 

We have analysed power-law models with power spectral index n — 2, 1, 0, —1 in 2D, with a cutoff in each case at the Nyquist 
wave number k c — 2048 fc/, where kf = 2n/Lhox is the fundamental mode associated with the box size. 

We choose ct(/cnl), the epoch when the scale 2-7r/fcNL is going nonlinear, as a measure of time. The first scale to go nonlinear 
is the one corresponding to the Nyquist wave number. This happens, by definition, when the variance a is unity. As a increases, 
successive larger scales enter in the nonlinear regime. The simulations were stopped at Anl = 2Z gr id, 4/ gr id, 8/ gr id, •■••, Lbox/2. 
In our study we have studied epochs till L b ox/16 goes nonlinear. 

The construction of h(x) uses measurements of the CPDF from several different epochs. We concentrated on smaller 
scales which are already in the highly nonlinear regime. Cell sizes of (4/ gr id) 2 , (8Z gr id) 2 and (16Z gr id) 2 were studied for four 
different epochs when Lbox/8, Lb ox /16, Lbox/32 and Lbox/64 have gone nonlinear. 

The growth rates of various modes in the linear regime were studied by Melott et al. (1988) for this PM code. The results 
at A = 3/ gr id arc equivalent to those obtained by a typical PM code at A = 8/ gr id, owing to the staggered mesh scheme. 
So we expect that our code performs well at the wavelength associated with four cells and since the collapse of 4Z gr id-size 
perturbations will give rise to condensations of diameter 2 Z gr id or less, the smallest cell size that can be safely resolved is 2 l ST id- 
The epochs we include in our study are not affected by the fact that simulations were started by Zel'dovich approximation. 
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Figure 3. The measured h(x) in 2D numerical simulations, as described in the text. The large dynamic range allows us to determine the 
function h(x) to an accuracy at the level of lO -3 , and order-of-magnitude increase compared to all previous studies in 2D or 3D. The 
function h(x) in general shows a power law profile for x < 1 which is followed by an exponential cutoff for x > 1. This exponential cutoff 
is more severe when there are more power at smaller scales. Since x = N/N c and N c is the typical occupation of cells inside clusters, this 
means that the CPDF Py(N) at any scale will have a power-law profile for N < N c and will have an exponential tail for N > 7V C . In an 
infinite catalogue this exponential tail will extend to infinitely large values of N. However, in a finite catalogue, Pi (TV) shows an abrupt 
cutoff at ./Vmax, due to the absence of any denser cells beyond that limit. Different length scales map different regions of h(x) curve. 
Previous studies have shown that effect of finite volume corrections increases with power on larger scales which reduces the available 
range of x for which h(x) can be measured accurately. The measured h(x) satisfies two integral constraints (Si = 52 = 1) and the power 
law index of x 2 h(x) should be exactly the same as the negative power-law index appearing in a(N c ). 

4.2 Data Analysis 

For analysing the data, computations of the count probability distribution function (CPDF) Pi(N) and 2CPDF Pi ir (Ni, N?) 
were performed by laying down a grid of mesh spacing I and counting the occupation number in each cell; this yields the 
probability of finding N particles in cell size I and also the joint probability of finding N\ and N2 particles in two different 
cells separated by a distance r. Statistics were improved by perturbing the grid in each orthogonal direction and keeping the 
mesh undistorted while repeating the counting process. We considered cells of size 4Z gr id, 8Z gr id and 16Z gr id in our studies. With 
these cell sizes we can reach probabilities as small as few times 10 -9 . We have also checked that our results do not change if 
we increase the number of sampling cells; we get exactly the same results even with ten times fewer cells. For computing the 
VPF we have only considered cell sizes where the condition Po > 1/e is satisfied; it was shown by Colombi et al. (1995) that 
the properties of voids larger than this are affected by the grid which is used to start simulations. 

To extend the range of N c for which we can study the scaling properties we have diluted the data while studying scaling 
properties of VPF. Different levels of dilution such as 2048 2 , 1024 2 , 512 2 and 256 2 were considered. Results of the VPF 
analysis are plotted in Figure 5, showing that the scaling ansatz predicts the VPF very well. We find a power-law profile for 
a(N c ) for large values of N c . At small values of N c , when the data is dominated by Poisson noise, a(N c ) tends to unity. Our 
results of analysis of VPF proves beyond doubt that the hierarchical ansatz is indeed a good approximation in the highly 
non-linear regime. Since a(N c ) acts as a generating function for the Sn parameters this means that all Sn parameters reach 
a constant value in this limit. Furthermore, as the Sn parameters are linear combinations of product of powers of amplitudes 
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Figure 4. The measured values of b(> x) from the 2D simulations described in the text. The bias b(> x) shows two different regimes 
similar to h(x): for small values of x, it increases slowly with x according to a power-law with index less than unity. This index is also 
related to the negative slope of cr(N c ) and slope of x 2 h(x) for x < 1. For larger values of x, b(> x) shows a steeper increase, proportional 
to x. For b(N) this means that moderately overdense cells show a slow increase of bias with cell occupancy N but for highly overdensc cells 
the bias increases linearly. The transition occurs for N = N c which is the typical occupancy of cells in over-dense regions. We have used 
two related but different methods to estimate bias. The first method is based on measuring cross-correlation of different classes of objects 
with respect to background mass distribution and the second method is based on computation of direct measurement of correlation 
function against the background mass distribution for different classes of objects. Good agreement between these two methods proves 
that the factorization property of bias and also shows that the bias for an overdense cell depends only on the intrinsic property of the 
object concerned and can be expressed as a unique function of x. 



of different vertices that appear in the tree-level correlation hierarchy, this will also mean that these vertices too are constant 
in highly non-linear regime. 

However, we have also found a significant departure from such scaling for n = —2 spectra. The maximum value of iV c 
for which the VPF can be probed reliably decreases with increasing power on larger scales and the simulation grid, which 
remains undistorted in underdense regions for spectra with more power at larger scales even in the highly non-linear regime, 
these two factors may be responsible for a spurious departure from scaling in this case. Strong motivation for suspecting that 
the apparent departure we detect in scaling for VPF for n = —2 may not be real is that we do not detect any corresponding 
effect in h(x). 

We use four neighbouring equidistant cells in a grid to evaluate the bias with respect to the central cell. We adopt two 
different methods to evaluate the bias which, as suggested by Bernardeau (1994), test the factorizability property of bias of 
over-dense cells predicted by the hierarchical ansatz. In the first method we compute the bias associated with each overdense 
cells with respect to background mass. This is done by casting the definition of bias in a following way: 



&(> N) = 



E 



N t >N 



P(Ni)b(N : 






1 /^E N1>N E N3 N 2 P(N U N 2 ) 



E Nl>N P(Ni) ) foV E Nl>N P(Ni) 

In the second method we directly compute the bias from correlation of two different over-dense cells: 



(79) 
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Figure 5. The measured values of cr(N c ) from 2D simulations described in the text. For a given cell size, the parameter N c increases as 
the system evolves and for a given epoch, the larger cells have higher values of N c . Cell size we analyze ranges from (4£ gr ; d ) 2 to (64Z gr ; d ) 2 . 
We have used different epochs of nonlinearity to construct the function <r(iV c ). We also use different levels of nonlinearity to increase 
available range of N c . Different dilution levels of the 4096 2 data, i.e. 2048 2 , 1024 2 , 512 2 and 256 2 were used to test the validity of the 
scaling ansatz. For power spectra n = 2, 1 and we do not detect any deviation from scaling models, but for n = — 2 there arc discernible 
departures. This could be due to various spurious effects which could influence the determination of cr(N c ), especially the initial grid 
which survives undistortcd in underdense regions if the initial spectrum is n = —2. Error estimation in determination of cr(N c ) was done 
by Colombi et al. (1995) who found that the error is proportional to 54 which increases with more power on larger scales. However all 
our data points satisfy the criteria Po > 1/e for avoiding errors introduced by the computational grid. This restricts the maximum cell 
size for a given epoch for which ff(N c ) can be computed. The other possibility for violation of particular scaling model which we have 
considered here arises from the fact that instead of Sjv parameters being constant at highly nonlinear regime they might show a very 
slow but steady increase with level of nonlinearity as reported by Lucchin et. al (1994) and Colombi ct. al (1996). However note that 
such an argument will also mean a clear departure of CPDF from scaling which we have not detected in our studies. 
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P(JVi)6(JVi 
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N t >N 



P(Ni 




(80) 



In both cases we compute the bias for cells whose occupancy is greater than some particular threshold x and then we vary 
the threshold to study 6(> x) for different values of x. We initially attempted to compute directly the function b(x) but we 
found measurements of the cumulative distribution 6(> x) to be much more stable. We also found that the bias computed 
from both methods mentioned above matches, and also the slope of the bias function for overdense cells agrees well with 
theoretical predictions. In particular, the bias increases with a; as a power law with the index of powerdaw being less than 
unity, whereas for highly overdense cells the bias is directly proportional to x. Determination of the bias is more difficult for 
spectra with more power on smaller scales, due to absence of long range correlations in such models. 

Computation of the 2VPF was done using a similar technique as in the computation of bias. We used four neighbouring 
cells for every cell to find the joint probability of finding the pair of cells to be devoid of any particles; we did this for several 
different separations of two cells. For a given value of the scaling parameter N c , we have also studied the scaling properties 
associated with the 2VPF for larger separation. Data was diluted to from 4096 2 to 2048 2 , 1024 2 , 512 2 and 256 2 to increase 
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Figure 6. We plot the quantity r 2 (N c )/N 2 described in the text for different cell separations r. We have shown results for only two 
different initial conditions for one particular cell size of (4/ gr id) 2 - For each cell size we consider different levels of dilutions. The lowest 
curve in each panel correspond to 4096 2 data. Each subsequent curve from bottom to top correspond to a dilution by a factor of f /4 in 
N, while the top-most curve correspond to 64 2 data. The solid lines at larger values of r correspond to values of t 2 (N c )/N 2 taken for 
scaling studies; this quantity exhibits scaling properties when plotted as a function of N c as shown in Figure 7. Departure from linear 
behaviour for small separation is due to terms higher order in (,ij/(,i which are neglected in our study of the 2VPF. 



Table 1. Parameters of the fitting function h(x) and b(x) in 2D 
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the available dynamic range of the scaling parameter N c as we did for the case of VPF. For large cell separations, where all 
higher order terms in £12 /£ are negligible, 
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(81) 



We found that the scaling function t(N c ) 2 /N 2 , which encodes the scaling properties of the 2VPF and also acts as a generating 
function for cumulant correlator (Munshi & Melott 1998) is described very accurately by the theoretical prediction, made on 
the basis of the hierarchical ansatz; see Figures 6 & 7. 

5 FINITE VOLUME CORRECTION 

The finite size of galaxy catalogues and numerical simulations always affects the determination of many-body correlation 
functions from real or simulated samples. The higher the order of correlation measured, the greater the contribution to it from 
exceptionally dense and consequently rare clumps. Their determination is therefore always affected by the presence or absence 
of such clumps in a finite sample. This manifests itself in a sharp cutoff in all the one-point functions we have defined such as 
P(N), b(N), and vi (N). Using the scaling arguments we have already described n it is however possible to supplement the 
measured quantities with the functions by h(x), b(x), and V2(x). These quantities have been computed by series expansion of 
multi-point VPF in increasing power of (£ij/£i) for second and third order by Bernardeau & Schaeffer (1992) and for fourth 
and fifth order in our Paper I. Using their results we show here that scaling arguments once verified at lower order in (£ij/£i) 
can be used to estimate finite volume errors for in estimates of the multi-point cumulant correlators. 

If Af max is the maximum occupancy of cells of certain length scales in a catalog, measured P{N) for N < iV max shows an 
abrupt cutoff for TV — iV max . We supplement the measured P(N) by n(^) for N > iV max to construct the corrected P C (N) 
(Munshi et al. 1998c); II(^) can be expressed in terms of the scaling function h(x) describing P(N): 



P C (N) 
P C {N) 



,-g-^w, 



11(1/) 



N < iV max 
iV > iV max . 



(82) 



The corrected CPDF satisfies two constraint equations similar to the case of the uncorrected CPDF. The constraints 
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Figure 7. The measured values of t 2 (N c )/N 2 from 2D simulations described in the text. We compute the 2VPF P(0, 0) and VPF -P(O) 
using a regular grid for different separation as explained in previous figures and use this to estimate t 2 (N c )/N 2 for different values of 
N c in the limit of large separations as described in the text. As in the case of VPF we repeat the process for different values of cell size 
from (4Zg r id) 2 to (64Z gr id) 2 - We also use different levels of dilution to increase the range of N c studied. We study 2VPF for 4096 2 , 2048 2 , 
1024 2 , 512 2 , 256 2 , 128 2 and 64 2 data points. All length scales which do not satisfy the criteria Pq < 1/e were excluded from analysis to 
avoid lattice effects. We also study different levels of nonlinearity and cell sizes as in the analysis of VPF. Dashed lines correspond to 
the theoretical prediction r 2 /y 2 oc y~( 1 + u} \ where —uj is also the slope of one point VPF. 



Y]p(JV) = l, ^P c (iV) = l, Y^iVP(iV) =N, ^2nP c (N) = N 
can be used to determine the constants A and B where 
A — H Q = / h(x)dx 



B — Hi = / xh(x)dx 

and 

S p c = E p -( P +l)i/ 1 S- 1 + ^. 

Similarly, the bias associated with overdense cells b(N) shows an abrupt cutoff at N > Af max which we supplement by its 
continuous analogue b{v), which can be expressed in terms of its scaling function b(x): 



(83) 

(84) 
(85) 

(86) 



P c (N)b c (N) = (l-C-V(N -■^±)\p(N)b(N) N<N„ 

P c {N)b c (N) = U(u)b(u) 



(87) 



N>N U 



where we have denoted ^2 N 2 P(N)b(N) by (N 2 }t,. The unknown parameters C and T> can be determined by using the 
fact that both b(N)P(N) and corrected b c (N)P c (N) satisfy similar constraints. For the zeroth moment: J2 b(N)P(N) = 0, 
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J2 b c (N)P COTT (N) = and similarly for the first moment we have £ Nb(N)P(N) = N, £ Nb c (N)P c (N) = iV. The constraint 
on the zeroth-order moment can be used to determine T>, and the first-order constraint fixes the value of C: 

/■oo 

C = Si = / xh(x)b{x)dx (88) 

•^iVmax/iV e 
/■OO 

£> = B = h(x)b(x)dx (89) 

•^iVmax/iV c 

and 

E£ c = (l-.B 1 )E£-(p+l)BoE£ +1 + ^. (90) 

Where we have defined the p th order moment of b(x)h(x) by B p and 



B p = / a: p &(a;)/i(:r)d:r. (91) 

Extending these calculations to the case of MCPDFs is fairly straightforward. The constraints ^ P(N)i> s (N) — Nv s , 
^2NP c (N)u^(N) = Nv a and ^P(A/> s (iV) = h T, pC ( N )»s( N ) = l can be used to determine the values of X and y 
appearing in the renormalization of P(N)u a (N) that corrects for the finite volume effect: 

P*(N)vZ(N) = (l-X-y(N-^^-)^P(N)v n (N) iV<iV max (92) 

P c (N)v c n (N) = U{v)v n {v) N>N n 



• max , 



where y and X are the first and zeroth order moments of h(x)u n (x) respectively and {N 2 ) Vn represents y^ N 2 P(N)v n (N)\ 
these quantities can be expressed by 

1 f°° 
X — Zf™ — — / xh(x)v n {x)dx (93) 

Vn JN max /N c 
1 f°° 

y = Z^ = — h{x)v n {x)dx (94) 



Vn 



k/JV c 



and 

E^ c = (1 - Z?)Y? - Z "»(p + 1)S^ + ^, (95) 

where we have denoted the p th moment of v n (x)h(x) by 

1 /- 00 
Zp n — — / x v h(x)v n {x]dx. (96) 

Vn JN maI /N c 

If we also notice that we have to correct N c for the finite volume effect using N£ = N c (l — 6HiY,3 + H2)- Incorporating 
this correction we can finally write 

c ^-(p+l^S^ + ff 

^ - (1-6H 1 S 3 +F2)^- 1 ( j 

yfc c (l-BQS^-b + l^oS^ + ff 

(1 - 6tfiE 3 + H 2 )p~i ( ^> 

(1 - Zj* )E- - Z- (p + 1)E£ , + J- 






(l-6 J ffiE 3 + /f 2 )P- 1 



We have kept only terms which are leading order in l/£ in the expressions and hence they are valid only in the highly nonlinear 
regime. We have also ignored a term in {N 2 ) Vs /S,2VnN 2 since we expect {N 2 ) Vn /v n N 2 to be of order unity. The inclusion of 
terms representing loops of lower-order diagrams in the calculations will bring not only moments of b(x)h(x) and v„(x)h(x) 
but also moments of powers of these scaling functions which we have neglected here. 

Cumulants and cumulant correlators have been measured with presently available N-body and galaxy catalogues. Although 
several analysis for finite volume correction for cumulants are available in the literature (Colombi et al. 1992, 1994, 1995, 
1996, Munshi et al. 1998), no such analysis has been done for cumulant correlators. Our analysis based on the hierarchical 
ansatz provides an estimate for finite volume correction and also can be used to correct such effects. Implementation of our 
method for measurements of MCCs will be presented elsewhere. 

It was shown by Szapudi & Szalay (1997) that cumulant correlators can actually be used to separate amplitudes associated 
with different tree topologies at low order. However, at higher orders the number of equations becomes less than the number 
of independent tree topologies rendering the system indeterminate. This can only be cured if we increase the number of points 
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by using multi-point cumulant correlators. This in principle will allow us not only to determine the amplitudes associated 
with tree topologies of arbitrary order but also allow an additional consistency check because the number of variables is 
less than the number of equations they have to satisfy. Of course going from two-point cumulant correlators to multi-point 
cumulant correlators makes numerical computation more complicated, especially because of the shape dependence of multi- 
point cumulant correlators for more than two points. With the availability of simulations with much larger dynamic range a 
determination of the v„ parameters will become possible in the reasonably near future. 

The procedure which we have introduced here for finite volume corrections however requires a priori knowledge of the v n 
parameters which one needs to determine from the same sample. This calls for an iterative procedure which can be applied 
until the whole process becomes convergent. 



6 CONCLUSIONS 

In this and the previous paper (Paper I), we have investigated properties of the multi-point cumulant correlators - the natural 
generalisations of the one-point cumulants and the two-point cumulant correlators. We also developed a method which can 
be used to determine the amplitudes associated with tree diagrams of different topologies, without making any simplistic 
assumptions about them. This method can be used to test any models of gravitational clustering in the highly non-linear 
regime, such as the hierarchical scaling ansatz. 

We have developed a diagrammatic method to represent the decomposition of multi-point cumulant correlators in terms 
of factorial moments and factorial correlators. Using this approach, we have found that every higher-order contribution to 
(Cij/C*) introduces a new function which corresponds to a new vertex in the tree- level hierarchy, labeled by v> n (N). It was 
shown in Paper I that these functions display similar scaling relations as the functions P(N) and b(N). We have used these 
scaling relations for v n (N) to estimate the effect of the finite size of N-body simulations or galaxy catalogues in estimates of 
statistical descriptors, such as cumulants, cumulant correlators and multi-point cumulant correlators. 

We have tested various predictions of the hierarchical ansatz, relating to one-point and two-point quantities in the highly 
non-linear regime, using high-resolution 2D numerical simulations. We find that the scaling functions h(x) and a(N c ) can 
describe the behaviour of the CPDF and VPF in highly non-linear regime for all length scales, as the hierarchical assumption 
predicts. These functions were, however, found to depend on initial conditions even in the highly nonlinear regime. As expected 
from general scaling arguments, it was found that a(N c ) attains a asymptotic power-law iV"" for large values of N c . However, 
such models are not able to predict the values of u>; more detailed dynamical arguments would have to be developed to make 
such a prediction. Numerically, however, we have found that lu depends on the form of the initial power spectrum, increasing 
with the relative amount of small-scale power. The index to also enters in the description of x 2 h(x) which we found to have 
a power law profile with power law index u> and an exponential cutoff for large values x. This transition occurs around the 
typical cell occupancy in overdense cells where N = N c , i.e. x — 1. 

Extending such studies of one-point quantities for the first time to their two-point analogues such as 2CPDF and 2VPF, 
we were able to check how overdense cells are biased with respect to the underlying mass distribution. We find, in accordance 
with hierarchical predictions, that the bias is an intrinsic property of collapsed objects. One can associate scaling variable x 
with every collapsed object, which is a function of cell occupancy TV, radius of the cell, and the variance of matter distribution 
£2 on that length scale. We also found that the bias is factorizable and for two different objects with different values of scaling 
parameters x\ and £2 and with intrinsic bias b(x\) and 6(0:2), it is possible to express their bias with respect to each other as 
6(iEi)&(a:2). 

Elsewhere in the literature, extensions of the Press-Schechter theory have been used to compute the bias of collapsed 
objects (Mo & White 1996) and also their higher-order one-point moments such as Sn parameters (Mo et al. 1997). These are 
in agreement with the results we have obtained, especially those presented in Paper I. However, it remains unclear whether 
the generalisation of this formalism to two-point quantities (such as the CCs) will produce accurate predictions. Other works 
which have focussed on predicting bias and higher-order correlation functions are related to extensions of the Zel'dovich 
approximation (Lee & Shandarin 1997,1998, Catelan et al. 1997). It would be a very important test of the usefulness of these, 
and other approximation schemes for non-linear gravitational evolution to see if they can actually reproduce the multi-point 
statistics of collapsed objects as well as the simple hierarchical ansatz we have discussed here. 

All these theoretical predictions were made in the lowest order of (£ij/£i), but were found to describe simulation results 
even when two cells are not separated by large distances. 

Just as the VPF is related to the generating function of the Sn parameters, the 2VPF is related to the generating function 
of the cumulant correlators (CCs). In earlier numerical studies (Munshi & Melott 1998) we have shown that, in general, CCs 
are factorizable, so that C pq — C p iC q i, when two cells are separated by a large distance. Such a decomposition of CCs is 
related to the factorization of the 2VPF which we have established in this paper, thus confirming our earlier findings. The 
factorization property of bias is also a direct consequence of the scaling properties we study here, so our results on the subject 
of bias also vindicate our earlier findings. 

Our numerical study was done in 2D and we plan to extend such analysis in case of 3D in near future, incorporating 
studies of higher-order correlations of over-dense cells. 

It is interesting to note that although we have demonstrated that the hierarchical ansatz (62) is a very good approximation 
to the behaviour of higher-order correlation functions, stable clustering is known to be violated in 2D in highly non-linear 
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regime (Munshi et al. 1998a). It is often argued that stable clustering is a necessary ingredient in the hierarchical ansatz. This 
appears not to be the case. We shall discuss this issue in more detail in a future paper. 

Finite sample (or simulation) size known to introduce a sharp cut-off in the CPDF for iV = A/max corresponding to the 
density of the densest densest cell in the catalog. Extending methods developed earlier (Munshi et al. 1998c) for correcting 
finite volume effects in the CPDF Pi (N) by supplementing the missing information beyond TV = iV max from scaling function 
h(x) (which can then be used to compute the Sn parameters), we developed a similar method for MCPDF in this paper. At the 
level of 2CPDF e.g., one can use the scaling property of b(N) which is encoded in scaling function b(x) to correct or estimate 
finite volume effects in the extraction of CCs. We show that this technique can be generalised in a straightforward manner 
to the case of MCPDFs. At present there is no other general formalism to estimate the finite volume error for multi-point 
statistics. Our method provides the first step in this direction. 
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